library(R.matlab)
library(ggplot2)
library(RColorBrewer)
library(reshape2)
library(multcomp)
library(nlme)
library(lme4)
library(ISLR)
library(scales)
library(svglite)
library(ez)
library(Hmisc)
library(dplyr)


rm(list=ls())

## --- User-defined parameters ---

savePlots = 0
plotAll = 0
genPlots = c(1, # 1) Average pupil traces; significance on No-Int conditions
             1) # 2) Residual Int - NoInt plots (+ perm test significance)

bPal = brewer.pal(4,'Dark2')

# Separate color schemes for WM and INT
condColors_WM = c('#802b00', '#cc4400', '#003399', '#3377ff') # Two-tone (orange family for aud, blue for vis)

## --- Other setup ---

numPlots = length(genPlots)
if (plotAll == 1) {
  genPlots = rep(1,numPlots)
}

# Set paths based on current machine
compname = Sys.info()
compname = compname[[4]]
if (compname == "DESKTOP-BHR0AU7") {
  loadPath = 'E:/ANL/Experiments/RESULTS/VAST/DATAforR/WholeTrial/'
  imgSaveDir = 'E:/ANL/Experiments/RESULTS/VAST/FIGURES/ms/'
} else if (compname == "MSI") {
  loadPath = 'E:/ANL/RESULTS/VAST/DATAforR/WholeTrial/'
  imgSaveDir = 'E:/ANL/RESULTS/VAST/FIGURES/ms/'
} else if (compname == "UMN202302853") {
  loadPath = 'D:/ANL/RESULTS/VAST/DATAforR/WholeTrial/'
  imgSaveDir = 'D:/ANL/RESULTS/VAST/FIGURES/ms/'
} else if (compname == "UMN202302892") {
  loadPath = 'E:/ANL/RESULTS/VAST/DATAforR/WholeTrial/'
  imgSaveDir = 'E:/ANL/RESULTS/VAST/FIGURES/ms/'
} else {
  stop('set directories for this machine')
}


## --- Load and format the pupillometry data ---

# AVERAGE DATA
pupil_avg_temp = readMat(paste(loadPath, 'avgPupil_wholeT.mat', sep=""))
# NA-pad the traces so I can put them all in a df
# first pass -- find the ENC and RET window lengths in each condition
ENC_lens = RET_lens = rep(0, length(pupil_avg_temp$timebases))
for (i in 1:length(pupil_avg_temp$timebases)) {
  ENC_lens[i] = sum(pupil_avg_temp$timebases[[i]][[1]] < 0)
  RET_lens[i] = sum(pupil_avg_temp$timebases[[i]][[1]] >= 0)
}
longestENC = max(ENC_lens)
longestRET = max(RET_lens)
# second pass -- nan pad
avg_loc = seq(2, 36, by=3) # indices of average pupil data
for (i in 1:length(pupil_avg_temp$timebases)) {
  pad_st = longestENC - ENC_lens[i]
  pad_end = longestRET - RET_lens[i]
  currpupil = pupil_avg_temp$avg.pupil[[avg_loc[i]]][[1]]
  currSEM = pupil_avg_temp$avg.pupil[[avg_loc[i]+1]][[1]]
  currpupil = c(matrix(data=NA, nrow=1, ncol=pad_st), currpupil)
  currpupil = c(currpupil, matrix(data=NA, nrow=1, ncol=pad_end))
  currSEM = c(matrix(data=NA, nrow=1, ncol=pad_st), currSEM)
  currSEM = c(currSEM, matrix(data=NA, nrow=1, ncol=pad_end))
  pupil_avg_temp$avg.pupil[[avg_loc[i]]][[1]] = currpupil
  pupil_avg_temp$avg.pupil[[avg_loc[i]+1]][[1]] = currSEM
}
# reconstruct a timebase
tStep = abs(pupil_avg_temp$timebases[[1]][[1]][2] - pupil_avg_temp$timebases[[1]][[1]][1])
TB_pupil = seq(0, (longestENC + longestRET - 1) * tStep, by = tStep)
TB_pupil = TB_pupil - TB_pupil[longestENC] # offset so start RET = t0

# Construct data frame
pupil_df = data.frame('time'=TB_pupil, 'AT_none'=pupil_avg_temp$avg.pupil[[11]][[1]], 'AT_none_SEM'=pupil_avg_temp$avg.pupil[[12]][[1]],
                      'AT_AT'=pupil_avg_temp$avg.pupil[[17]][[1]], 'AT_AT_SEM'=pupil_avg_temp$avg.pupil[[18]][[1]],
                      'AT_AS'=pupil_avg_temp$avg.pupil[[14]][[1]], 'AT_AS_SEM'=pupil_avg_temp$avg.pupil[[15]][[1]],
                      'AS_none'=pupil_avg_temp$avg.pupil[[2]][[1]], 'AS_none_SEM'=pupil_avg_temp$avg.pupil[[3]][[1]],
                      'AS_AT'=pupil_avg_temp$avg.pupil[[8]][[1]], 'AS_AT_SEM'=pupil_avg_temp$avg.pupil[[9]][[1]],
                      'AS_AS'=pupil_avg_temp$avg.pupil[[5]][[1]], 'AS_AS_SEM'=pupil_avg_temp$avg.pupil[[6]][[1]],
                      'VT_none'=pupil_avg_temp$avg.pupil[[29]][[1]], 'VT_none_SEM'=pupil_avg_temp$avg.pupil[[30]][[1]],
                      'VT_AT'=pupil_avg_temp$avg.pupil[[35]][[1]], 'VT_AT_SEM'=pupil_avg_temp$avg.pupil[[36]][[1]],
                      'VT_AS'=pupil_avg_temp$avg.pupil[[32]][[1]], 'VT_AS_SEM'=pupil_avg_temp$avg.pupil[[33]][[1]],
                      'VS_none'=pupil_avg_temp$avg.pupil[[20]][[1]], 'VS_none_SEM'=pupil_avg_temp$avg.pupil[[21]][[1]],
                      'VS_AT'=pupil_avg_temp$avg.pupil[[26]][[1]], 'VS_AT_SEM'=pupil_avg_temp$avg.pupil[[27]][[1]],
                      'VS_AS'=pupil_avg_temp$avg.pupil[[23]][[1]], 'VS_AS_SEM'=pupil_avg_temp$avg.pupil[[24]][[1]])

# AVERAGE RESIDUALS
ENC_lens = RET_lens = rep(0, length(pupil_avg_temp$timebases.res))
for (i in 1:length(pupil_avg_temp$timebases.res)) {
  ENC_lens[i] = sum(pupil_avg_temp$timebases.res[[i]][[1]] < 0)
  RET_lens[i] = sum(pupil_avg_temp$timebases.res[[i]][[1]] >= 0)
}
longestENC = max(ENC_lens)
longestRET = max(RET_lens)
# second pass -- nan pad
avg_loc = seq(2, 24, by=3) # indices of average pupil data
for (i in 1:length(pupil_avg_temp$timebases.res)) {
  pad_st = longestENC - ENC_lens[i]
  pad_end = longestRET - RET_lens[i]
  currpupil = pupil_avg_temp$avg.pupil.res[[avg_loc[i]]][[1]]
  currSEM = pupil_avg_temp$avg.pupil.res[[avg_loc[i]+1]][[1]]
  currpupil = c(matrix(data=NA, nrow=1, ncol=pad_st), currpupil)
  currpupil = c(currpupil, matrix(data=NA, nrow=1, ncol=pad_end))
  currSEM = c(matrix(data=NA, nrow=1, ncol=pad_st), currSEM)
  currSEM = c(currSEM, matrix(data=NA, nrow=1, ncol=pad_end))
  pupil_avg_temp$avg.pupil.res[[avg_loc[i]]][[1]] = currpupil
  pupil_avg_temp$avg.pupil.res[[avg_loc[i]+1]][[1]] = currSEM
}
# reconstruct a timebase
tStep = abs(pupil_avg_temp$timebases.res[[1]][[1]][2] - pupil_avg_temp$timebases[[1]][[1]][1])
TB_pupil_res = seq(0, (longestENC + longestRET - 1) * tStep, by = tStep)
TB_pupil_res = TB_pupil_res - TB_pupil_res[longestENC] # offset so start RET = t0

# Construct data frame
pupil_res_df = data.frame('time'=TB_pupil_res, 
                          'AT_AT'=pupil_avg_temp$avg.pupil.res[[11]][[1]], 'AT_AT_SEM'=pupil_avg_temp$avg.pupil.res[[12]][[1]],
                          'AT_AS'=pupil_avg_temp$avg.pupil.res[[8]][[1]], 'AT_AS_SEM'=pupil_avg_temp$avg.pupil.res[[9]][[1]],
                          'AS_AT'=pupil_avg_temp$avg.pupil.res[[5]][[1]], 'AS_AT_SEM'=pupil_avg_temp$avg.pupil.res[[6]][[1]],
                          'AS_AS'=pupil_avg_temp$avg.pupil.res[[2]][[1]], 'AS_AS_SEM'=pupil_avg_temp$avg.pupil.res[[3]][[1]],
                          'VT_AT'=pupil_avg_temp$avg.pupil.res[[23]][[1]], 'VT_AT_SEM'=pupil_avg_temp$avg.pupil.res[[24]][[1]],
                          'VT_AS'=pupil_avg_temp$avg.pupil.res[[20]][[1]], 'VT_AS_SEM'=pupil_avg_temp$avg.pupil.res[[21]][[1]],
                          'VS_AT'=pupil_avg_temp$avg.pupil.res[[17]][[1]], 'VS_AT_SEM'=pupil_avg_temp$avg.pupil.res[[18]][[1]],
                          'VS_AS'=pupil_avg_temp$avg.pupil.res[[14]][[1]], 'VS_AS_SEM'=pupil_avg_temp$avg.pupil.res[[15]][[1]])

rm(pupil_avg_temp)

# Load in permutation testing results
pupil_stat_temp = readMat(paste(loadPath, 'pupilStats_wholeT.mat', sep=""))
stat_tb = unlist(pupil_stat_temp$perm.TBs[1])

# Approach: Code time by "y", where y will be come the vertical position of the
# significance bars. If the permutation test was not significant, y will be NA
# at that time point. For significant results, the actual value of y will be
# offset to fit appropriately on the plot axes.

stat_results = data.frame()
# AT Intervening task
stats_temp = pupil_stat_temp$perm.results[[2]]
for (ind in seq(1,11,by=2)) {
  comparison = unlist(stats_temp[ind])
  signif_samples = stats_temp[[ind+1]][[1]]
  y = rep(NA, length(stat_tb))
  for (sr in 1:nrow(signif_samples)) {
    y[signif_samples[sr,1]:signif_samples[sr,2]] = 0
  }
  stat_results = rbind(stat_results,
                       data.frame('time' = stat_tb, 'y' = y) %>%
                         mutate(int_task = "AT",
                                comparison = comparison))
}
# Add the AS Intervening task
stats_temp = pupil_stat_temp$perm.results[[3]]
for (ind in seq(1,11,by=2)) {
  comparison = unlist(stats_temp[ind])
  signif_samples = stats_temp[[ind+1]][[1]]
  y = rep(NA, length(stat_tb))
  if (nrow(signif_samples) > 0) {
    for (sr in 1:nrow(signif_samples)) {
      y[signif_samples[sr,1]:signif_samples[sr,2]] = 0
    }
  }
  stat_results = rbind(stat_results,
                       data.frame(time = stat_tb, y = y) %>%
                         mutate(int_task = "AS",
                                comparison = comparison))
}

# Cleanup
# only keep time points >0 (WM retention / Intervening task interval)
stat_results = stat_results %>% dplyr::filter(time >= 0)

# increment y values so the lines are above the pupil traces and don't overlap
# with one another
# ...AT intervening task
uniq_comps_AT = unique(stat_results$comparison[stat_results$int_task=="AT"])
for (c in 1:length(uniq_comps_AT)) {
  offset = (c-1) * 0.1
  stat_results$y[!is.na(stat_results$y) & 
                 stat_results$comparison == uniq_comps_AT[c]] = 1.8 - offset
}
# ...AS intervening task
uniq_comps_AS = unique(stat_results$comparison[stat_results$int_task=="AS"])
for (c in 1:length(uniq_comps_AS)) {
  offset = (c-1) * 0.1
  stat_results$y[!is.na(stat_results$y) & 
                   stat_results$comparison == uniq_comps_AS[c]] = 1.8 - offset
}

# Make a separate DF containing the comparison labels
# Note: one comparison omitted (VT WM vs. VS WM) because it didn't reach
# significance for either the AT or AS Intervening task.
stat_labels = data.frame(int_task = c(rep("AT", 5), rep("AS", 5)),
                         comparison = c(uniq_comps_AT[1:5],
                                        uniq_comps_AS[1:5]),
                         time = rep(-1,10),
                         y = rep(seq(1.8,1.4,by=-0.1), 2))


# *************************************************************************************
# _____________________________________________________________________________________
#
# --- PLOTTING ---  
# _____________________________________________________________________________________
# *************************************************************************************

##  1) Average pupil traces; significance on No-Int conditions
if (genPlots[1]==1) {
  
  # --- No INT TASK ---
  tplot = ggplot(data=pupil_df, aes(x=time)) +
    geom_segment(x=0, y=-Inf, xend=0, yend=Inf, size=2, linetype=2, color='#999999ff') +
    geom_segment(x=2, y=-Inf, xend=2, yend=Inf, size=2, linetype=2, color='#999999ff') +
    geom_segment(x=-Inf, y=0, xend=Inf, yend=0, size=2, linetype=2, color='#999999ff') +
    
    # ribbon plots for SEM
    geom_ribbon(aes(ymin=VS_none - VS_none_SEM, ymax=VS_none + VS_none_SEM), fill=condColors_WM[4], alpha=0.2) +
    geom_ribbon(aes(ymin=VT_none - VT_none_SEM, ymax=VT_none + VT_none_SEM), fill=condColors_WM[3], alpha=0.2) +
    geom_ribbon(aes(ymin=AS_none - AS_none_SEM, ymax=AS_none + AS_none_SEM), fill=condColors_WM[2], alpha=0.2) +
    geom_ribbon(aes(ymin=AT_none - AT_none_SEM, ymax=AT_none + AT_none_SEM), fill=condColors_WM[1], alpha=0.2) +
    
    # add on group plots
    geom_line(aes(y=VS_none), color=condColors_WM[4], size=3) + 
    geom_line(aes(y=VT_none), color=condColors_WM[3], size=3) +
    geom_line(aes(y=AS_none), color=condColors_WM[2], size=3) +
    geom_line(aes(y=AT_none), color=condColors_WM[1], size=3) +
    
    scale_x_continuous(name='Time (s)', breaks=seq(-2, 5, by=1), expand=c(0,0)) +
    scale_y_continuous(name='Pupil Diameter (z)', breaks=seq(-0.4, 1, by=0.2), expand=c(0,0)) +
    coord_cartesian(xlim=c(-2,5.5), ylim=c(-0.55,1.15)) +
    theme(panel.background = element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.major.x = element_blank(),
          aspect.ratio = 1,
          panel.border = element_rect(colour = "black", fill = NA, size=3.5),
          plot.margin=unit(c(2,2,2,2), "lines"),
          legend.direction = 'vertical',
          axis.text.x = element_text(size=40, colour = "black"),
          axis.text.y = element_text(size=40, colour = "black"),
          axis.title.x = element_blank(),
          axis.title.y = element_blank())
  
  windows()
  print(tplot)
  
  if (savePlots == 1) {
    ggsave(filename = paste(imgSaveDir,'wholeT_pupil_noINT.svg',sep=""), width=10, height=10, units="in")
    graphics.off()
  }
  
  # --- AT INT TASK ---
  tplot = ggplot(data=pupil_df, aes(x=time)) +
    geom_segment(x=0, y=-Inf, xend=0, yend=Inf, size=2, linetype=2, color='#999999ff') +
    geom_segment(x=2, y=-Inf, xend=2, yend=Inf, size=2, linetype=2, color='#999999ff') +
    geom_segment(x=-Inf, y=0, xend=Inf, yend=0, size=2, linetype=2, color='#999999ff') +
    
    # ribbon plots for SEM
    geom_ribbon(aes(ymin=VS_AT - VS_AT_SEM, ymax=VS_AT + VS_AT_SEM), fill=condColors_WM[4], alpha=0.2) +
    geom_ribbon(aes(ymin=VT_AT - VT_AT_SEM, ymax=VT_AT + VT_AT_SEM), fill=condColors_WM[3], alpha=0.2) +
    geom_ribbon(aes(ymin=AS_AT - AS_AT_SEM, ymax=AS_AT + AS_AT_SEM), fill=condColors_WM[2], alpha=0.2) +
    geom_ribbon(aes(ymin=AT_AT - AT_AT_SEM, ymax=AT_AT + AT_AT_SEM), fill=condColors_WM[1], alpha=0.2) +
    
    # add on group plots
    geom_line(aes(y=VS_AT), color=condColors_WM[4], size=3) + 
    geom_line(aes(y=VT_AT), color=condColors_WM[3], size=3) +
    geom_line(aes(y=AS_AT), color=condColors_WM[2], size=3) +
    geom_line(aes(y=AT_AT), color=condColors_WM[1], size=3) +
    
    scale_x_continuous(name='Time (s)', breaks=seq(-2, 5, by=1), expand=c(0,0)) +
    scale_y_continuous(name='Pupil Diameter (z)', breaks=seq(-0.4, 1, by=0.2), expand=c(0,0)) +
    coord_cartesian(xlim=c(-2,5.5), ylim=c(-0.55,1.15)) +
    theme(panel.background = element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.major.x = element_blank(),
          aspect.ratio = 1,
          panel.border = element_rect(colour = "black", fill = NA, size=3.5),
          plot.margin=unit(c(2,2,2,2), "lines"),
          legend.direction = 'vertical',
          axis.text.x = element_text(size=40, colour = "black"),
          axis.text.y = element_text(size=40, colour = "black"),
          axis.title.x = element_blank(),
          axis.title.y = element_blank())
  
  windows()
  print(tplot)
  
  if (savePlots == 1) {
    ggsave(filename = paste(imgSaveDir,'wholeT_pupil_atINT.svg',sep=""), width=10, height=10, units="in")
    graphics.off()
  }
  
  # --- AS INT TASK ---
  tplot = ggplot(data=pupil_df, aes(x=time)) +
    geom_segment(x=0, y=-Inf, xend=0, yend=Inf, size=2, linetype=2, color='#999999ff') +
    geom_segment(x=2, y=-Inf, xend=2, yend=Inf, size=2, linetype=2, color='#999999ff') +
    geom_segment(x=-Inf, y=0, xend=Inf, yend=0, size=2, linetype=2, color='#999999ff') +
    
    # ribbon plots for SEM
    geom_ribbon(aes(ymin=VS_AS - VS_AS_SEM, ymax=VS_AS + VS_AS_SEM), fill=condColors_WM[4], alpha=0.2) +
    geom_ribbon(aes(ymin=VT_AS - VT_AS_SEM, ymax=VT_AS + VT_AS_SEM), fill=condColors_WM[3], alpha=0.2) +
    geom_ribbon(aes(ymin=AS_AS - AS_AS_SEM, ymax=AS_AS + AS_AS_SEM), fill=condColors_WM[2], alpha=0.2) +
    geom_ribbon(aes(ymin=AT_AS - AT_AS_SEM, ymax=AT_AS + AT_AS_SEM), fill=condColors_WM[1], alpha=0.2) +
    
    # add on group plots
    geom_line(aes(y=VS_AS), color=condColors_WM[4], size=3) + 
    geom_line(aes(y=VT_AS), color=condColors_WM[3], size=3) +
    geom_line(aes(y=AS_AS), color=condColors_WM[2], size=3) +
    geom_line(aes(y=AT_AS), color=condColors_WM[1], size=3) +
    
    scale_x_continuous(name='Time (s)', breaks=seq(-2, 5, by=1), expand=c(0,0)) +
    scale_y_continuous(name='Pupil Diameter (z)', breaks=seq(-0.4, 1, by=0.2), expand=c(0,0)) +
    coord_cartesian(xlim=c(-2,5.5), ylim=c(-0.55,1.15)) +
    theme(panel.background = element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.major.x = element_blank(),
          aspect.ratio = 1,
          panel.border = element_rect(colour = "black", fill = NA, size=3.5),
          plot.margin=unit(c(2,2,2,2), "lines"),
          legend.direction = 'vertical',
          axis.text.x = element_text(size=40, colour = "black"),
          axis.text.y = element_text(size=40, colour = "black"),
          axis.title.x = element_blank(),
          axis.title.y = element_blank())
  
  windows()
  print(tplot)
  
  if (savePlots == 1) {
    ggsave(filename = paste(imgSaveDir,'wholeT_pupil_asINT.svg',sep=""), width=10, height=10, units="in")
    graphics.off()
  }
}


##  2) Residual Int - NoInt plots (+ perm test significance)
if (genPlots[2]==1) {
  
  # --- AT INT TASK ---
  tplot = ggplot(data=pupil_res_df, aes(x=time)) +
    geom_segment(x=0, y=-Inf, xend=0, yend=Inf, size=2, linetype=2, color='#999999ff') +
    geom_segment(x=2, y=-Inf, xend=2, yend=Inf, size=2, linetype=2, color='#999999ff') +
    geom_segment(x=-Inf, y=0, xend=Inf, yend=0, size=2, linetype=2, color='#999999ff') +
    
    # ribbon plots for SEM
    geom_ribbon(aes(ymin=VS_AT - VS_AT_SEM, ymax=VS_AT + VS_AT_SEM), fill=condColors_WM[4], alpha=0.2) +
    geom_ribbon(aes(ymin=VT_AT - VT_AT_SEM, ymax=VT_AT + VT_AT_SEM), fill=condColors_WM[3], alpha=0.2) +
    geom_ribbon(aes(ymin=AS_AT - AS_AT_SEM, ymax=AS_AT + AS_AT_SEM), fill=condColors_WM[2], alpha=0.2) +
    geom_ribbon(aes(ymin=AT_AT - AT_AT_SEM, ymax=AT_AT + AT_AT_SEM), fill=condColors_WM[1], alpha=0.2) +
    
    # add on group plots
    geom_line(aes(y=VS_AT), color=condColors_WM[4], size=3) + 
    geom_line(aes(y=VT_AT), color=condColors_WM[3], size=3) +
    geom_line(aes(y=AS_AT), color=condColors_WM[2], size=3) +
    geom_line(aes(y=AT_AT), color=condColors_WM[1], size=3) +
    
    # add on significance labels and bars
    geom_line(aes(x=time, y=y, group=comparison), 
              data=stat_results %>% dplyr::filter(int_task == "AT")) +
    
    # add on labels for significance bars
    geom_label(aes(x=time, y=y, label=comparison),
               data = stat_labels %>% dplyr::filter(int_task=="AT"), 
               size=4) +
    
    scale_x_continuous(name='Time (s)', breaks=seq(-2, 5, by=1), expand=c(0,0)) +
    scale_y_continuous(name='Pupil Diameter Difference(z)', breaks=seq(-0.4, 1, by=0.2), expand=c(0,0)) +
    # coord_cartesian(xlim=c(-2,5.5), ylim=c(-0.55,1.15)) +
    coord_cartesian(xlim=c(-2,5.5), ylim=c(-0.55,1.9)) +
    theme(panel.background = element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.major.x = element_blank(),
          aspect.ratio = 1,
          panel.border = element_rect(colour = "black", fill = NA, size=3.5),
          plot.margin=unit(c(2,2,2,2), "lines"),
          legend.direction = 'vertical',
          axis.text.x = element_text(size=40, colour = "black"),
          axis.text.y = element_text(size=40, colour = "black"),
          axis.title.x = element_blank(),
          axis.title.y = element_blank())
  
  windows()
  print(tplot)
  
  if (savePlots == 1) {
    ggsave(filename = paste(imgSaveDir,'wholeT_pupil_res_atINT.svg',sep=""), width=10, height=10, units="in")
    graphics.off()
  }
  
  # --- AS INT TASK ---
  tplot = ggplot(data=pupil_res_df, aes(x=time)) +
    geom_segment(x=0, y=-Inf, xend=0, yend=Inf, size=2, linetype=2, color='#999999ff') +
    geom_segment(x=2, y=-Inf, xend=2, yend=Inf, size=2, linetype=2, color='#999999ff') +
    geom_segment(x=-Inf, y=0, xend=Inf, yend=0, size=2, linetype=2, color='#999999ff') +
    
    # ribbon plots for SEM
    geom_ribbon(aes(ymin=VS_AS - VS_AS_SEM, ymax=VS_AS + VS_AS_SEM), fill=condColors_WM[4], alpha=0.2) +
    geom_ribbon(aes(ymin=VT_AS - VT_AS_SEM, ymax=VT_AS + VT_AS_SEM), fill=condColors_WM[3], alpha=0.2) +
    geom_ribbon(aes(ymin=AS_AS - AS_AS_SEM, ymax=AS_AS + AS_AS_SEM), fill=condColors_WM[2], alpha=0.2) +
    geom_ribbon(aes(ymin=AT_AS - AT_AS_SEM, ymax=AT_AS + AT_AS_SEM), fill=condColors_WM[1], alpha=0.2) +
    
    # add on group plots
    geom_line(aes(y=VS_AS), color=condColors_WM[4], size=3) + 
    geom_line(aes(y=VT_AS), color=condColors_WM[3], size=3) +
    geom_line(aes(y=AS_AS), color=condColors_WM[2], size=3) +
    geom_line(aes(y=AT_AS), color=condColors_WM[1], size=3) +
    
    # add on significance labels and bars
    geom_line(aes(x=time, y=y, group=comparison), 
              data=stat_results %>% dplyr::filter(int_task == "AS")) +
    
    # add on labels for significance bars
    geom_label(aes(x=time, y=y, label=comparison),
               data = stat_labels %>% dplyr::filter(int_task=="AS"), 
               size=4) +
    
    scale_x_continuous(name='Time (s)', breaks=seq(-2, 5, by=1), expand=c(0,0)) +
    scale_y_continuous(name='Pupil Diameter Difference(z)', breaks=seq(-0.4, 1, by=0.2), expand=c(0,0)) +
    # coord_cartesian(xlim=c(-2,5.5), ylim=c(-0.55,1.15)) +
    coord_cartesian(xlim=c(-2,5.5), ylim=c(-0.55,1.9)) +
    theme(panel.background = element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.major.x = element_blank(),
          aspect.ratio = 1,
          panel.border = element_rect(colour = "black", fill = NA, size=3.5),
          plot.margin=unit(c(2,2,2,2), "lines"),
          legend.direction = 'vertical',
          axis.text.x = element_text(size=40, colour = "black"),
          axis.text.y = element_text(size=40, colour = "black"),
          axis.title.x = element_blank(),
          axis.title.y = element_blank())
  
  windows()
  print(tplot)
  
  if (savePlots == 1) {
    ggsave(filename = paste(imgSaveDir,'wholeT_pupil_res_asINT.svg',sep=""), width=10, height=10, units="in")
    graphics.off()
  }
}



